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ABSTRACT 

Motivation: Cellular information processing can be described 
mathematically using differential equations. Often, external 
stimulation of cells by compounds such as drugs or hormones 
leading to activation has to be considered. Mathematically, the 
stimulus is represented by a time-dependent input function. 

Parameters such as rate constants of the molecular interactions 
are often unknown and need to be estimated from experimental data, 
e.g. by maximum likelihood estimation. For this purpose, the input 
function has to be defined for all times of the integration interval. 
This is usually achieved by approximating the input by interpolation 
or smoothing of the measured data. This procedure is suboptimal 
since the input uncertainties are not considered in the estimation 
process which often leads to overoptimistic confidence intervals of 
the inferred parameters and the model dynamics. 
Results: This article presents a new approach which includes 
the input estimation into the estimation process of the dynamical 
model parameters by minimizing an objective function containing 
all parameters simultaneously. We applied this comprehensive 
approach to an illustrative model with simulated data and compared it 
to alternative methods. Statistical analyses revealed that our method 
improves the prediction of the model dynamics and the confidence 
intervals leading to a proper coverage of the confidence intervals of 
the dynamic parameters. The method was applied to the JAK-STAT 
signaling pathway. 

Availability: MATLAB code is available on the authors' website 
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1 INTRODUCTION 

Mathematical modeling of biological systems has become a widely 
used approach to better understand the system behavio r as a whole 
rather than observing isolated parts (iKitanol I2002I). The rapid 
devel opment in quantitative molecular biology ( Cox and Mannl 
l201ll) enables to calibrate mathematical models to experimental 
data and therefore generating model predictions. For parametric 
models, such as ordinary differential equation (ODE) models, 
calibration can be performed e.g. by maximum likelihood estimation 
(MLE). Uncertainties of the measurements generated by biological 
variability or technical limitations have to be considered in the 
calibration process as they propagate to the parameter estimates. 



An accurate method for calculat ing confidence int ervals is given by 
the profile likelihood approach feaue et a/.ll2009h . 

One unresolved issue in data-based modeling is the insufficient 
consideration of input measurement uncertainties in the parameter 
estimation process. In this article, we solve this problem with a 
new approach. The method presented in the following includes 
the input estimation into the estimation process of the dynamical 
model parameters by minimizing an objective function containing 
all parameters simultaneously. An illustrative model with known 
parameter values and noise distributions reveals that our method is 
able to (i) increase the precision of the parameter estimates and (ii) 
correct the coverage of the likelihood-based confidence intervals. 

Furthermore, we applied the method to a model of th e JAK- 
STAT signaling pathway introduced by ISwameve et al ] d2003h 
where quantitative biological data are available, leading to different 
estimated dynamics of the model trajectories compared to the 
standard method. 



2 METHODS 

2.1 Models of biochemical reaction networks 

A mathematical model of a biochemical reaction network is given by a 
set of coupled ODEs 

x(t)=f(*(t)M),p) (1) 



y(t)=g(x(t),u(t),s) + €(t) With 6;(0=A/"(0,or2), 



(2) 



where xt(t) are the internal states of the system, m(t) are the time dependent 
external stimuli, also called input functions and yi(t) are the noisy observables 
of the system with a noise level oi . 

The set of parameters that defines M. uniquely contains dynamic 
parameters p, i.e. for instance rate constants, scaling and offset parameters 
s as well as the initial concentrations x(to) 

9 = \p,x(t 0 ),s}eR + . (3) 

All analyses have been performed on a logarithmic scale. 

2.2 Parameter estimation 

We performed MLE which is equivalent to the weighted least squares method 
for Gaussian noise. Thus, an objective function 

2 



N I 
i=l \ 



' y] -y{x(f u 0)Mti\0) 



is minimized and 



0 = argminx (0). 



(4) 



(5) 
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is the maximum likelihood estimator where y] denotes the measured data 
and y(x(ti,0),u(ti),0) represents the model output. In order not to rely on 



© The Author(s) 2012. Published by Oxford University Press. 

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.0rg/licenses/by/3.O), which 
permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 



M.Schelker et al. 



the mere point estimate, likelihood profiles will be computed as discussed in 
more detail in the Supplementary Information. 



input parametrization even though the spline is smoothed by penalizing the 
curvature. 



2.3 Existing methods 

For integration of the ODEs, a continuous representation of a discretely 
measured input is required. In the past, time-depe ndent input functions u(t ) 
have been represented by constant i nput functions iBachmann et all 1201 1|) ; 
linear interpolations of the dataset lOuach et all 120071 : ISwameve et all 
1 2003b , cubi c interpolation splines or by cubic smoothing splines ( Maiwald 
and Timmer. 1200a iRaue et all\200$) . ~~~~ 

If the functional representation of the input data is estimated separately 
from the dynamical model parameters, the deviations from the true 
underlying dynamics of the input are not considered in the estimation of 
the model parameters. This can result on the one hand in biased estimates 
of the dynamic parameters and on the other hand, even more likely, in 
underestimated confidence intervals as the input uncertainties have not been 
taken into account. This problem shall be addressed in the following by 
defining a suitable input parametrization that enables to estimate input and 
dynamical parameters comprehensively. 

2.4 Splines 

A spline is a piecewise polynomial that interpolates a given set of data points 
y\,...,y^ N at time points t\,...,t^ in a smooth and continuous manner. The 
v-values to be interpolated are called control points and the corresponding t- 
values knots. Boundary conditions determine either the slope or the curvature 
of the spline at the first and the last data point. Here, cubic splines with natural 
boundary conditions, i.e. the curvature at both end points is set to zero, are 
utilized. 

A cubic interpolation spline can be described either by a set of knots and 
the corresponding four coefficients for each polynomial piece or by a linear 
combination of a set of basis functions in the basis spline (B -spline) notation. 
For a knot sequence t\ , . . . ,fry and i ts corr esponding control points v\ , . . . , v#, 
the B-spline is given by dde Booilll978l) 



N-p-l 



S(t)= ^2 v i b i,p(t) with te[ti +p ,t N - p ], 



(6) 



where N is the number of data points and p is the degree of the basis functions 
bi^ p it), thus, in case of cubic polynomials, p = 3. The closed-form expression 
of the basis functions is given in the Supplementary Information. 

Another type of spline, that is frequently used, is the cubic smoothing 
spline. A smoothing spline 



<S(r) = argmin 



N 



I d 2 S(t) 



dt 2 



dt 



(7) 



can, depending on its respective measurement uncertainties 07, deviate 
from the data points yj in order to form a smoother curve. Therefore, 
an objective function is to be minimized, containing both the deviation 
from the data and the curvature, i.e. the second term in J7J. The so- 
called smoothing or penalization parameter X e [0, 00) controls the trade-off 
between interpolating the noisy data and obtaining a smooth function. As 
a consequence, X = 0 results in an interpolation spline, whereas X —> 00 
corresponds to a linear fit of the input data. 

Interpolation splines often lead to strongly oscillating input function 
estimates. Therefore, smoothing splines usually lead to more realistic 
outcomes if the smoothing parameter X can be chosen reasonably. 

A smoothing spline can be unambiguously represented either by the data 
points y [ , . . . , y^ , the knot sequence t\,...,t^ and the smoothing parameter X 
or by the knots t\ , . . . , f# and the corresponding control points vi , . . . , v#, that 
are the obtained by evaluation of |Equation 7| i.e. v z - = «S(f/). Therefore, for the 
approach presented in the following, the control points v\ , . . . , v# are used as 



2.5 Assessment criteria 

The performance of the approach is assessed by the accuracy 

Ai:=\(p0i))\ with ft=log 10 (£/ftrue,/). (8) 
and the precision 

Vi'.=Jv*r<j>0i% (9) 
i.e. the standard deviation of the distribution of the parameter estimates, if 
many independent estimations are performed over simulated data for known 
true parameter values. 

Moreover, the deviation from the true trajectory is quantified by the score 

N 2 

S(yi0 = ^(^^^ue)-^^)) ' ( 10 ) 
7=1 

where tj e [^o^end] is densely sampled at 500 equidistant time points. 

The coverage ratio, i.e. the probability that the true parameters are inside 
the confidence intervals, is calculated using the likelihood-based confidence 
intervals for many noise realization as described in more detail in the 
Supplementary Information. This value can be compared to the desired 
coverage, i.e. the corresponding a level. 

2.6 Approach 

The new approach aims to include the estimation of the input function into 
the maximum likelihood estimation of the model parameters. Therefore, a 
parametrization for the input measurements is needed. 

A cubic spline is uniquely defined by the set of data points it is fitted to. 
This is best visible in the B-spline notation introduced in |Equation 6| Thus, 
the input can be parametrized using the control points v\ , . . . , of the spline 
as parameters. 

The optimal set of parameters is then obtained by minimizing an objective 
function 

2 



x\e*y 



^^( yl-yi{x(tj,e%u(tj,e%e^ 
/=u=i\ Gij 



* {d 2 u k (t,e*)\ 



(id 



dt 



that contains both input and downstream parameters as well as a term that 
penalizes the curvature of the input functions. The smoothing parameter X is 
selected based on x 2 -statistic as described in the Supplementary Information. 
Note that the first term in |Equation 1 1 | accounts for the deviation of the model 
ytitj) from the data yi for all observables including the input. 
The extended parameter set to be optimized is 



6* = {p,x(t 0 ),v,s}. 



(12) 



In the following, 6*^6 will be used to keep the notation as simple as 
possible. 

Input measurements are usually concentrations or relative measurements 
like receptor activation and therefore strictly positive. A spline can become 
negative between the positive measurement points leading to reversed 
reaction fluxes and incorrect model behavior. To avoid negativity of the 
spline function, the input is transformed to logarithmic scale. 



3 APPLICATIONS 

In the following, the novel approach, which will also be called 
comprehensive approach, is compared to the standard approach, 
where the input is estimated from the measurement data before 
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Fig. 1. Schematic representation of the illustrative model. The forward 
reaction from Xi to X2 is catalyzed by the input u. This reaction is reversible 
with a rate constant k-\. Furthermore, X2 reacts to X3 with rate constants 
&2 and k-2 

the model parameters are calibrated. In lSection 3. II the smoothing 
parameter X was set to zero for both the standard and the 
comprehensive approach. By contrast, for both approaches in 
ISection 3?2\ the smoothing parameter was chosen optimally to 
the measured data based on x 2_s tatistics as described in the 
Supplementary Information. 

3.1 An illustrative model 

The first application of the new approach is a toy model with three 
dynamic states Xi , X2 and X3 and one input variable u. As depicted 
in |FigureT| the model consists of two reversible reactions where the 
first forward reaction (Xi — > X2) is catalyzed by a time-dependent 
input function u(t). The second forward reaction (X2 — > X3) as well 
as both reverse reactions (X3 — > X2 and X2 — > Xi ) are considered 
to be first-order mass action kinetics and do not depend on the input 
directly. 

The ODE system describing the dynamics of the model is given 
in the Supplementary Information. The observables are chosen 
to be yi = logio(xi) for 1 = 1,2,3 and J4 = logio(w) as for protein 
measurement techniques log-normally distr ibuted measurement 
uncertainties are assumed (iKreutz et alWlwh . 

For the simulation, N = 12 time points between t$ = 0 and ? enc j = 
50 were chosen with an equidistant sampling. The kinetic rate 
constants were set to k\=0.01, k_ \ = 1 , &2 = 0.5 and k_2 = 0.\. 
The initial concentrations were chosen close to the steady states, i.e. 
x\ (0) = 30, *2(0) — 20 and x 3 (0) = 50, as this maximizes the impact 
of the input on the variables. For the underlying true input dynamics, 
a Gaussian function 

9 10 3 if'-25\ 2 
u(t) = 10~ 2 + — — e~ 2 (—) (13) 
4V2tt 

was chosen which has its maximum w(7 max )~ 100 at t mSLX = 25. 

In order to evaluate the assessment criteria, M = 2000 noise 
realizations were simulated with a noise level ax = 0.1 for the 
compounds X/ and <r u = 03 for the input u(t). A comparison 
of several different noise levels is given in the Supplementary 
Information. For each noise realization, the parameters were 
estimated as described in Section algorithm in the Supplementary 
Information. 

In Panel A of |Figure~2| the model trajectories are plotted for 
one noise realization. The identifiability of the model parameters 
was analyzed by the profile likelihood approach as described in 
the Supplementary Information. The likelihood-based confidence 
intervals were calculated for multiple confidence levels 

ae{0.05,0.1,... ,0.95, 0.99} (14) 

to investigate the coverage. 



3.1.1 Evaluation of the assessment criteria The performance 
of the novel approach was challenged by the three assessment 
criteria introduced in ISection 231 In Panel B of |Figure 2| the 
parameter estimate distributions are shown. In terms of accuracy, 
both approaches perform similarly well as the dashed lines indicating 
the mean values are not distinguishable. The precision, by contrast, 
is increased by a factor of 1.04 up to 1.48 when the novel method 
is applied. This can also be seen in the plot where all parameter 
distributions in blue are narrower than the green ones. The values 
of A and V as well as the ratio of these values for both approaches 
are given in the Supplementary Information. 

Panel C of |Figure 2| shows that for all four observables, the 
deviation from the estimated to the true trajectory is significantly 
smaller in case of the comprehensive approach. 

In Panel D of |Figure 2[ the coverage ratio is plotted versus 
the confidence level for both approaches. In case of the standard 
approach, the measurement uncertainties of the input are not taken 
into account. As a consequence, the confidence intervals, at least for 
some of the parameters, are too small. This results in a coverage 
ratio that is smaller than the corresponding confidence interval 
as indicated by the deviation of the green solid line from the 
black dashed line in Panel D of |Figure~2] Underestimating the 
size of a confidence interval can lead to false-positive rejections of 
specific values for the respective parameter and therefore represents 
a serious issue. For the comprehensive approach by contrast, the 
coverage ratio almost equals the confidence level for all seven 
model parameters. According to this, the uncertainty of the input 
measurements is appropriately translated into parameter confidence 
intervals. 

Apparently, the first two rate constants k\ and k_ \ are affected the 
most by the choice of the estimation approach. This results from the 
model structure as the fluxes of the first reaction X^ =^X2 depend 
directly on the input. Consistently, the coverage of the rate constants 
&2 and k_2 of the second reaction X2 ^ X3 are barely influenced 
by the input and its uncertainties. 

3.2 The JAK-STAT signaling pathway 

For the second application of the novel approach, a data-based 
model of JAK2-STAT5 signaling pathway dSwameve et all 120031) 
was investigated. In this model the downstream system depends 
on a measured input. The model structure is depicted in |Figure 3| 
The ODEs describing the dynamics of the model are given in the 
Supplementary Inf o rmatio n. 

Swameye et al. d2QQ3h performed quantitative immunoblotting 
in order to investigate the Epo-induced activation of the JAK2- 
STAT5 signaling pathway. Due to experimental limitations, only 
sums of several species were measurable leading to the following 
observables: 

y 1 = s 1 • (*2 + 2x 3 ) , yi = $2 • (* 1 + *2 + 2x 3 ) . 

y 3 =s 3 -u, 

where y\ is the total amount of phosphorylated STAT5 in 
the cytoplasm, y^ the total amount i.e. unphosphorylated and 
phosphorylated STAT5 in the cytoplasm and y 3 is the measurement 
of the input function i.e. the amount of the Epo-induced 
phosphorylation of EpoR. Furthermore, the immunoblotting method 
could only provide relative protein concentrations. Therefore, 



i531 



M.Schelker et al. 




1.5 2 2.5 3 3.5 3 3.5 4 4.5 5 5.5 6 6.5 7 0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1 



logto(s) confidence level a 



Fig. 2. (A) Logarithmic plot of model trajectories and simulated data. (B) Probability densities for the estimates of the model parameters. (C) Score distributions. 
(D) The coverage ratio is plotted against the confidence level to check the size of the confidence intervals. The dashed line indicates identity. In all figures, the 
green solid line represents the trajectory of the standard approach and the blue one of the comprehensive approach. The black lines indicate the true dynamics 



scaling parameters s\ to 53 are required and the observables are 
given in arbitrary units (au). 

In order to calibrate the model to experimental data, a set of 
initial conditions for all parameters is needed. First, the initial 
concentrations of all phosphorylated model species were set to zero: 

[pSTAT5](0) -> 0, [pSTAT5 2 ](0) -> 0, 

[npSTAT5 2 ](0)^0. 

Furthermore, in order to resolve structural identifiability issues, 
the initial value of STAT5 and the scaling parameter of the input 
were set to one: 

[STAT5](0)->1, 53^1. 

Thereby, the forward rate constant k\ is expressed in units of the 
input pEpoR and the total amount of STAT5, which is conserved in 
this system, is set to one. In a second step, a set of initial guesses 
for the remaining parameters was selected using Latin hypercube 
sampling (LHS). Subsequently, the parameters were estimated by 
maximum likelihood estimation. 

Panel A of [FT gure 4 1 shows the resulting model trajectories and 
the corresponding time course data. The model variables, which 
correspond to the optimal set of parameters for each approach, are 
depicted in |Figure 5] For the model species STAT5, npSTAT52 and 
for the input pEpoR, the trajectories differ only sightly. However, 
for pSTAT5 and pSTAT52, a completely different scale is observed 
for both approaches. This results from parameter values of ki and 
&3 that differ strongly which will be discussed later on. 

The identifiability of the model parameters was analyzed by 
the profile likelihood approach as described in the Supplementary 
Information. The resulting likelihood profiles are shown in Panel B 
of |Figure"4] The corresponding values of the parameter estimates 
and their confidence intervals are given in lTable ll 




Fig. 3. Schematic representation of the JAK-STAT signaling pathway. 
The phosphorylated Erythropoietin receptor (pEpoR) acts as input for the 
signaling cascade and induces phosphorylation of STAT5. Phosphorylated 
STAT5 (pSTAT5) then dimerizes (pSTAT5 2 ) and translocates to the nucleus 
(npSTAT5 2 ). There it can act as transcription factor and thereby induce 
protein synthesis. Furthermore, with some delay r, the pSTAT5 dimers 
dissociate to unphosphorylated ST AT5 monomers and tran slocate back to 
the cytoplasm. (Figure adapted from lSwameve et al] 120031) ) 

The parameter estimates of t he rate constants Jc\ and £4 can be 
directly compared to results of iRaue et al\ d20Q9h . However, the 
point estimates obtained by the standard approach do not exactly 
match the previously obta ined values . This is due to the choice of 
the input representation: in lRaue etal. 1 d2009h . a smoothing spline of 
the input data was used without considering the input uncertainties 
for the weighting of the data points and for the determination 
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time [min] logi 0 (#) 

Fig. 4. (A) Plot of model trajectories and experimental data. (B) Likelihood profiles for the model parameters. The red dashed line indicates the confidence 
threshold for a = 0.95. The green solid line represents the likelihood profile calculated with the standard approach and the blue one with the comprehensive 
approach 
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Fig. 5. The green solid line represents the trajectory calculated with the 
standard approach and the blue one with the comprehensive approach. 
Note the different scaling of the left (standard approach) and the right 
(comprehensive approach) j-axis for pSTAT5 and pSTAT52 

of the smoothing parameter. In this article by contrast, the input 
uncertainties have been estimated from several datasets as described 
in the Supplementary Information. This allows for a weighting of the 
data according to the first term in Equation (7) and for a data-b ased 
selection of the smoothing parameter X. In iRaue et ali d2009h . the 
initial concentration of STAT5 (x\(0)) and the scaling parameters 
s\ and S2 were structurally non-identifiable. In our article, the 
parameter value of [STAT5](0) was set to one leading to identifiable 
scaling parameters having different physical units. Therefore, these 
parameter values are not comparable with previous results. The 
results for the rate constants &2 and &3 are discussed later on in 
this article. 

3.2.1 Comparison of the results Here, the results of the standard 
approach and the comprehensive approach shall be discussed 
and compared. This analysis requires a different strategy than in 
ISection 3TT1 as the underlying truth is not know for this system. 

At first, it can be observed that the qualitative behavior of the fits 
in Panel A of |Figure "4| differs between both approaches. The input 
observable 373, for instance, has a small plateau at about 25 min 
in case of the comprehensively estimated input. This is due to 
downstream information, as in the total concentration of pSTAT5, 
i.e. the observable y\, this plateau is also present with some delay. 
Furthermore, for the standard approach, the model trajectories of y\ 
and deviate significantly from the data at the last measurement 
time point, i.e. £ = 60 min. This is due to the input measurement 



which has a higher value at the last time point (t = 60 min) than for 
the proceeding measurement time points (t = 40, 50 min). In contrast, 
in case of the comprehensive approach, the experimental data 
including the uncertainties of all three observables were considered, 
leading to model trajectories that better reflect the measured data. 

One, at first glance, confusing outcome of the model calibration is 
the difference in the parameter estimates for &2 and £3 as depicted by 
the likelihood profiles in Panel B of |Figure"4| These two parameters 
show a completely opposed behavior: for the standard approach, 
the process of dimerization, which is described by ki, is faster than 
estimated by the comprehensive approach. In contrast, £3, which 
describes the translocation of the pSTAT5 dimers to the nucleus, is 
smaller for the standard approach than the estimate resulting from 
the comprehensive approach. As a consequence, the trajectories 
of the dynamic variables pSTAT5 and pSTAT52 show a different 
scaling for both estimation approaches. Therefore, in |Figure'5| the 
concentrations obtained by the standard approach are indicated by 
the y-axis on the left and those of the comprehensive approach by 
the right j-axis. 

T his opposing behavior of fa and fa was also observed bv Raue 
etal. (120091) . who investigated the variability of the trajectories when 
plotting all parameter sets along the likelihood profile that are within 
the point- wise 95% confidence interval. This results from the fact 
that only the sum of pSTAT5 and pSTAT52 is measured but not their 
individual contribution. 

Here, the switching from faster dimerization and slower 
translocation to the opposed configuration can be explained in 
terms of the / 2 -landscape. For the comprehensive approach, a 
local minimum can be identified by LHS, which corresponds to 
the situation of the standard approach as depicted in Panels A and 
B of |Figure~6| Although both minima are in sufficient agreement 
with the data because x 2 differs by much less than the number of 
fitted parameters, the optimal point estimates change. The parameter 
values of the minima are given in the Supplementary Information. 
However, by optimizing the input in the comprehensive approach, 
the x 2 -landscape changes and another minimum with the switched 
parameter configuration becomes optimal. 

The selection of a different local minimum in case of 
the comprehensive approach also leads to the counterintuitive 
observation that the confidence intervals for the parameters £4, 
5i and S2 are larger for the standard approach than for the novel 
method. 
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LHS for the comprehensive approach PL(&2) PL(^3) 




sorted fits logio(#) logio(^) 

Fig. 6. (A) Analysis of local minima for the comprehensive approach. Several local minima, indicated by dashed lines, have been identified by LHS. B and C: 
Likelihood profiles for ki and £3 . The green solid line represents the likelihood profile calculated with the standard approach and the blue one represents the 
profile around the optimum for the comprehensive approach. The red dashed line and the red asterisks denote a selected local minimum for the comprehensive 
approach with the switched parameter configuration 

Table 1. Parameter estimates and likelihood-based 95% confidence intervals for both approaches 



Parameter [units] Standard approach Comprehensive approach 







6 


confib 


conf ub 


6 


confib 


conf ub 


k\ [min -1 


nM -1 -5 -1 ] 


2.79 


2.18 


3.70 


1.95 


1.27 


3.33 


&2 [min -1 


nM -1 ] 


1.00 xlO 5 


4.04 


00 


0.11 


0.08 


0.14 


k 3 [min -1 ] 




0.12 


0.09 


0.17 


9.84 xlO 4 


0.34 


00 


&4 [min -1 




2.19 


1.69 


3.12 


1.49 


1.20 


2.01 


5i [nM -1 ] 




1.34 


1.24 


1.46 


1.25 


1.15 


1.36 


s 2 [nM -1 ] 




1.00 


0.95 


1.05 


0.95 


0.91 


0.99 



Hence, the result of the comprehensive approach reveals that an 
unambiguous estimation of both rates is not feasible based on the 
measured data. Therefore, additional experiments would be required 
in order to clarify which parameter configuration corresponds to the 
real biological situation. 



4 CONCLUSIONS 

It has been shown that the presented approach for estimating model 
parameters and input comprehensively provides several advantages. 
For a toy model with simulated data, the three assessment criteria 
defined in I Section 2. 5 1 have been evaluated for different noise level 
combinations of input and downstream observables. The analysis 
revealed that the novel approach outperforms the standard method in 
terms of (i) precision, (ii) deviation from the truth and (iii) coverage 
of the confidence intervals. The accuracy of the parameter estimates 
does not differ significantly for both approaches. The performance 
is best when the measurement uncertainties of the downstream 
observables are smaller than those of the input measurements. 

By applying the new method to a model of the JAK-STAT 
signaling pathway, it could be shown that also for realistic 
experimental studies, one benefits from the comprehensive input 
estimation. The optimal model trajectories indicate that downstream 
information can be used in order to improve the fit of the input and 
vice versa. 
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